Primary Analyses
# start by pulling in imputed data
Load(imp_prim_aim2)
# first need to make clinic.female a factor in all iterations
long <- complete(imp_prim_aim2, action = "long", include = TRUE)
long$clinic.female <- factor(long$clinic.female)
# now make it back into a mids object
imp_prim_aim2 <- as.mids(long)
# create function for pooling mm function
# first need to create tidy dataframe with 6 colums: term, estimate, std.error, statistic, p.value, and nobs (degrees of freedom) from MMLogit class (from GEE version)
# rows is the number of rows in the combined dataset
fitlist.mm <- function(object, rows) {
# first get summary
sum <- summary(object)
# now make tidy version
sum_tidy <- as_tibble(sum$mean.table) %>%
add_column(term = rownames(sum$mean.table), .before = 1) %>%
mutate(nobs = rows) %>%
rename(estimate = Estimate,
std.error = `Model SE`,
statistic = `Chi Square`,
p.value = `Pr(>Chi)`)
sum_tidy
}
# now make glanced list with rows n.clusters, max.cluster.size, alpha, and nobs
glanced.mm <- function(object, rows) {
data.frame(n.clusters = unname(object$control["n_subj"]),
max.cluster.size = unname(object$control["max_n_visit"]),
alpha = unname(object$alpha),
nobs = rows)
}
# barnard rubin function pulled directly from mice code https://rdrr.io/cran/mice/src/R/barnard.rubin.R
barnard.rubin <- function(m, b, t, dfcom = 999999) {
lambda <- (1 + 1 / m) * b / t
lambda[lambda < 1e-04] <- 1e-04
dfold <- (m - 1) / lambda^2
dfobs <- (dfcom + 1) / (dfcom + 3) * dfcom * (1 - lambda)
dfold * dfobs / (dfold + dfobs)
}
pool.mm <- function(object, rows) {
call <- paste0("pool(object = ")
# based on pool.fitlist from https://rdrr.io/cran/mice/src/R/pool.R
p <- map_df(object, ~ fitlist.mm(.x, rows))
# we don't want group_by to change the order of the terms
p <- p %>% mutate(term = factor(term, levels = unique(term)))
pooled <- p %>%
group_by(term) %>%
summarise(
m = n(),
qbar = mean(estimate),
ubar = mean(std.error^2),
b = var(estimate),
t = ubar + (1 + 1 / m) * b,
dfcom = rows,
df = barnard.rubin(m, b, t, dfcom),
riv = (1 + 1 / m) * b / ubar,
lambda = (1 + 1 / m) * b / t,
fmi = (riv + 2 / (df + 3)) / (riv + 1)
)
pooled <- data.frame(pooled)
names(pooled)[names(pooled) == "qbar"] <- "estimate"
# based on https://rdrr.io/cran/mice/src/R/get.df.R
g <- map_df(object, ~ glanced.mm(.x, rows))
# finally need make mipo object to work with mice summary function
mipo.mm <- list(m = pooled$m[1],
pooled = pooled,
glanced = g)
class(mipo.mm) <- c("mipo", "data.frame")
mipo.mm
}
# write updated GEE with MD correction function using geeglm instead of gee
# adpated from https://rdrr.io/cran/geesmv/src/R/GEE.var.md.R
geeglm.var.md <- function(formula, family, data, corstr){ # could not get id to work without writing into the function, so not transferable...
# get model design matrix
m <- model.frame(formula, data)
mt <- attr(m, "terms")
data$response <- model.response(m, "numeric")
mat <- as.data.frame(model.matrix(formula, m))
# fit GEE model
gee.fit <- geeglm(formula = formula,
family = family,
id = clinic.female, ### UPDATE ID HERE IF MOVING TO DIFFERENT MODEL WITH DIFFERENT ID ###
data = data,
corstr = corstr)
beta_est <- gee.fit$coefficient
alpha <- gee.fit$geese$alpha # updated this line to get alpha directly
len <- length(beta_est)
len_vec <- len^2
### Estimate the robust variance for \hat{\beta}
data$id_fit <- unname(gee.fit$id) # updated with unname and updated variable to "id_fit"
cluster<-cluster.size(data$id_fit)
ncluster<-max(cluster$n)
size<-cluster$m
mat$subj <- rep(unique(data$id_fit), cluster$n)
if(is.character(corstr)){
var <- switch(corstr,
"independence"=cormax.ind(ncluster),
"exchangeable"=cormax.exch(ncluster, alpha),
"AR-M"=cormax.ar1(ncluster, alpha),
"unstructured"=summary(gee.fit)$working.correlation,)
}else{
print(corstr)
stop("'working correlation structure' not recognized")
}
if(is.character(family)){
family <- switch(family,
"gaussian"="gaussian",
"binomial"="binomial",
"poisson"="poisson")
}else{
if(is.function(family)){
family <- family()[[1]]
}else{
print(family)
stop("'family' not recognized")
}
}
# calculate MD correlation
cov.beta<-unstr<-matrix(0,nrow=len,ncol=len)
step11<-matrix(0, nrow=len, ncol=len)
for (i in 1:size){
y<-as.matrix(data$response[data$id_fit==unique(data$id_fit)[i]])
covariate<-as.matrix(subset(mat[,-length(mat[1,])], mat$subj==unique(data$id_fit)[i]))
var_i=var[1:cluster$n[i],1:cluster$n[i]]
if (family=="gaussian"){
xx<-t(covariate)%*%solve(var_i)%*%covariate
step11<-step11+xx
}else if (family=="poisson") {
D<-mat.prod(covariate, exp(covariate%*%beta_est))
Vi <- diag(sqrt(c(exp(covariate%*%beta_est))),cluster$n[i])%*%var_i%*%diag(sqrt(c(exp(covariate%*%beta_est))),cluster$n[i])
xx<-t(D)%*%solve(Vi)%*%D
step11<-step11+xx
}else if (family=="binomial"){
D<-mat.prod(covariate, exp(covariate%*%beta_est)/((1+exp(covariate%*%beta_est))^2))
Vi <- diag(sqrt(c(exp(covariate%*%beta_est)/(1+exp(covariate%*%beta_est))^2)),cluster$n[i])%*%var_i%*%diag(sqrt(c(exp(covariate%*%beta_est)/(1+exp(covariate%*%beta_est))^2)),cluster$n[i])
xx<-t(D)%*%solve(Vi)%*%D
step11<-step11+xx
}
}
step12<-matrix(0,nrow=len,ncol=len)
step13<-matrix(0,nrow=len_vec,ncol=1)
step14<-matrix(0,nrow=len_vec,ncol=len_vec)
p<-matrix(0,nrow=len_vec,ncol=size)
for (i in 1:size){
y<-as.matrix(data$response[data$id_fit==unique(data$id)[i]])
covariate<-as.matrix(subset(mat[,-length(mat[1,])], mat$subj==unique(data$id_fit)[i]))
var_i=var[1:cluster$n[i],1:cluster$n[i]]
if (family=="gaussian"){
xy<-t(covariate)%*%solve(var_i)%*%solve(cormax.ind(cluster$n[i])-covariate%*%solve(step11)%*%t(covariate)%*%solve(var_i))%*%(y-covariate%*%beta_est)
step12<-step12+xy%*%t(xy)
step13<-step13+vec(xy%*%t(xy))
p[,i]<-vec(xy%*%t(xy))
}else if (family=="poisson") {
D<-mat.prod(covariate, exp(covariate%*%beta_est))
Vi <- diag(sqrt(c(exp(covariate%*%beta_est))),cluster$n[i])%*%var_i%*%diag(sqrt(c(exp(covariate%*%beta_est))),cluster$n[i])
xy<-t(D)%*%solve(Vi)%*%solve(cormax.ind(cluster$n[i])-D%*%solve(step11)%*%t(D)%*%solve(Vi))%*%(y-exp(covariate%*%beta_est))
step12<-step12+xy%*%t(xy)
step13<-step13+vec(xy%*%t(xy))
p[,i]<-vec(xy%*%t(xy))
}else if (family=="binomial"){
D<-mat.prod(covariate, exp(covariate%*%beta_est)/((1+exp(covariate%*%beta_est))^2))
Vi <- diag(sqrt(c(exp(covariate%*%beta_est)/(1+exp(covariate%*%beta_est))^2)),cluster$n[i])%*%var_i%*%diag(sqrt(c(exp(covariate%*%beta_est)/(1+exp(covariate%*%beta_est))^2)),cluster$n[i])
xy<-t(D)%*%solve(Vi)%*%solve(cormax.ind(cluster$n[i])-D%*%solve(step11)%*%t(D)%*%solve(Vi))%*%(y-exp(covariate%*%beta_est)/(1+exp(covariate%*%beta_est)))
step12<-step12+xy%*%t(xy)
step13<-step13+vec(xy%*%t(xy))
p[,i]<-vec(xy%*%t(xy))
}
}
for (i in 1:size){
dif<-(p[,i]-step13/size)%*%t(p[,i]-step13/size)
step14<-step14+dif
}
cov.beta<-solve(step11)%*%(step12)%*%solve(step11)
cov.var<-size/(size-1)*kronecker(solve(step11), solve(step11))%*%step14%*%kronecker(solve(step11), solve(step11))
return(list(cov.beta=diag(cov.beta), cov.var=cov.var))
}
# now create my formulas # FIGURE OUT KEY OUTCOME...
ua_form <- formula(Y_imp ~ group.female)
# given different outcomes
a_form <- formula(Y_imp ~ group.female + rcs(day.female, 3) + rcs(age.female, 3) + district.female + who_stage
+ edu_cat.female + edu_cat.male + rcs(soc_sup_ps.female, 3) + rcs(soc_sup_ns.female, 3))
cc_form <- formula(Y_raw ~ group.female + rcs(day.female, 3) + rcs(age.female, 3) + district.female + who_stage
+ edu_cat.female + edu_cat.male + rcs(soc_sup_ps.female, 3) + rcs(soc_sup_ns.female, 3))
lc_form <- formula(Y_lc ~ group.female + rcs(day.female, 3) + rcs(age.female, 3) + who_stage
+ edu_cat.female + edu_cat.male + rcs(soc_sup_ps.female, 3) + rcs(soc_sup_ns.female, 3))
wc_form <- formula(Y_wc ~ group.female + rcs(day.female, 3) + rcs(age.female, 3) + who_stage
+ edu_cat.female + edu_cat.male + rcs(soc_sup_ps.female, 3) + rcs(soc_sup_ns.female, 3))Unadjusted Analysis
# run MM across all iterations of the imputation
prim_ua <- complete(imp_prim_aim2, "all") %>%
map(~ mm(mean.formula = ua_form,
lv.formula = ~ 1,
id = clinic.female,
q = 120,
data = .x))
# nw use my own pooling and summary functions to create output value
prim1 <- prim_ua %>%
pool.mm(., nrow(complete(imp_prim_aim2, 1))) %>%
summary
prim1 term estimate std.error statistic df p.value
1 (Intercept) -0.3739615 0.2853536 -1.3105196 919.2391 0.1903475
2 group.femaleIntervention -0.2244781 0.4155680 -0.5401719 894.9972 0.5892129# figure out sigma (this will influence what model I use moving forward)
map_dbl(prim_ua, ~.x$alpha) 1 2 3 4 5 6 7
0.20715753 0.19304287 0.16686178 0.10544682 0.19883016 0.13653610 0.13335233
8 9 10 11 12 13 14
0.18525253 0.12000055 0.19267406 0.20116210 0.20421226 0.20146474 0.22414458
15 16 17 18 19 20 21
0.18951001 0.18037471 0.22575955 0.21998425 0.14993807 0.20113387 0.14683238
22 23 24 25 26 27 28
0.25943309 0.14674950 0.18827007 0.26428377 0.18206040 0.11662212 0.23410446
29 30 31 32 33 34 35
0.21602462 0.16475053 0.17627552 0.14246570 0.17934063 0.08745458 0.19754740
36 37 38 39 40 41 42
0.20756357 0.22380212 0.22294152 0.15245710 0.15387600 0.19473212 0.12214312
43 44 45 46 47 48 49
0.20908635 0.16538781 0.16958884 0.14111852 0.22287817 0.20982342 0.19384049
50
0.18111764 Min. 1st Qu. Median Mean 3rd Qu. Max.
0.08745 0.15281 0.18889 0.18219 0.20746 0.26428 Adjusted Analysis
# run MM across all iterations of the imputation
prim_aa <- complete(imp_prim_aim2, "all") %>%
map(~ mm(mean.formula = a_form,
lv.formula = ~ 1,
id = clinic.female,
q = 120,
data = .x)) # now use my own pooling and summary functions to create output value
prim2 <- prim_aa %>%
pool.mm(., nrow(complete(imp_prim_aim2, 1))) %>%
summary
prim2 term estimate
1 (Intercept) -1.5925091461
2 group.femaleIntervention 0.2767165489
3 rcs(day.female, 3)day.female 0.0005555342
4 rcs(day.female, 3)day.female' -0.0003977941
5 rcs(age.female, 3)age.female 0.0261540237
6 rcs(age.female, 3)age.female' -0.0141429747
7 district.femaleInhassunge 1.4102873742
8 district.femaleNamacurra 0.6696444163
9 district.femaleMocubela 1.8155757242
10 district.femaleMaganja da Costa 1.0263414931
11 district.femaleGilé 1.9182074848
12 district.femaleQuelimane -1.6532366627
13 who_stageII 0.1447331499
14 who_stageIII 0.0381904072
15 who_stageIV -4.8916154140
16 edu_cat.femaleSome Primary School (Grades 1-7) 0.4241532740
17 edu_cat.femaleCompleted Primary School (Grade 7) 0.8350862827
18 edu_cat.femaleSome Secondary School (Grades 8-10) 0.7101020602
19 edu_cat.femaleCompleted Secondary School (Grade 10) 0.8008523636
20 edu_cat.femaleCollege/Higher Education 0.2940432026
21 edu_cat.maleSome Primary School (Grades 1-7) -0.3877111432
22 edu_cat.maleCompleted Primary School (Grade 7) -0.3238327695
23 edu_cat.maleSome Secondary School (Grades 8-10) -0.4125946159
24 edu_cat.maleCompleted Secondary School (Grade 10) -0.5640782535
25 edu_cat.maleCollege/Higher Education 0.0702739600
26 rcs(soc_sup_ps.female, 3)soc_sup_ps.female 0.0197512548
27 rcs(soc_sup_ps.female, 3)soc_sup_ps.female' -0.0168043446
28 rcs(soc_sup_ns.female, 3)soc_sup_ns.female -0.0515028820
29 rcs(soc_sup_ns.female, 3)soc_sup_ns.female' 0.0876116032
std.error statistic df p.value
1 1.515312e+00 -1.050944654 519.9711 0.29377218
2 4.679028e-01 0.591397507 895.7170 0.55440333
3 6.667471e-04 0.833200821 459.1421 0.40516469
4 9.108560e-04 -0.436725604 452.2277 0.66251880
5 4.639395e-02 0.563737861 268.7346 0.57340267
6 7.133156e-02 -0.198270935 261.4593 0.84298717
7 6.588707e-01 2.140461525 902.9085 0.03258481
8 6.114241e-01 1.095220822 919.2259 0.27370678
9 7.187809e-01 2.525909746 880.8233 0.01171422
10 8.968500e-01 1.144384835 581.9247 0.25293456
11 7.993255e-01 2.399782761 918.7612 0.01660300
12 1.592929e+00 -1.037859783 854.3141 0.29962900
13 2.851372e-01 0.507591178 390.3731 0.61202673
14 5.073564e-01 0.075273339 191.4814 0.94007582
15 7.144477e+02 -0.006846709 928.9133 0.99453863
16 2.360172e-01 1.797128437 367.4615 0.07313623
17 3.725015e-01 2.241833287 490.7698 0.02541816
18 3.409925e-01 2.082456473 659.8314 0.03768544
19 4.920298e-01 1.627649998 673.3683 0.10406690
20 6.823149e-01 0.430949383 274.6260 0.66684323
21 2.954143e-01 -1.312431833 542.5533 0.18992940
22 3.501676e-01 -0.924793771 610.0156 0.35543887
23 3.420418e-01 -1.206269608 531.5551 0.22825016
24 4.240538e-01 -1.330204429 414.7491 0.18418191
25 4.385138e-01 0.160254831 621.4103 0.87273242
26 4.863692e-02 0.406095945 346.8205 0.68492255
27 6.305428e-02 -0.266506006 337.6657 0.79001220
28 4.348321e-02 -1.184431371 279.3873 0.23724897
29 5.117247e-02 1.712084553 340.4323 0.08779176Alternative Adjusted Analysis
# run GEE across all iterations of the imputation (updated formula that works with md correction below)
prim_alt_aa <- complete(imp_prim_aim2, "all") %>%
map(~ geeglm(formula = Y_imp ~ group.female + rcs(day.female, 3) + rcs(age.female, 3)
+ edu_cat.female + edu_cat.male + rcs(soc_sup_ps.female, 3) + rcs(soc_sup_ns.female, 3),
family = binomial,
id = clinic.female,
data = .x,
corstr = "exchangeable"))
# can pool with mice functions for GEE
prim3 <- prim_alt_aa %>%
pool() %>%
summary()
prim3 term estimate
1 (Intercept) -0.4054243186
2 group.femaleIntervention -0.1943709051
3 rcs(day.female, 3)day.female 0.0006313232
4 rcs(day.female, 3)day.female' -0.0004617857
5 rcs(age.female, 3)age.female 0.0183105443
6 rcs(age.female, 3)age.female' 0.0011467792
7 edu_cat.femaleSome Primary School (Grades 1-7) 0.2153147648
8 edu_cat.femaleCompleted Primary School (Grade 7) 0.7094574826
9 edu_cat.femaleSome Secondary School (Grades 8-10) 0.6064644735
10 edu_cat.femaleCompleted Secondary School (Grade 10) 0.7853264552
11 edu_cat.femaleCollege/Higher Education 0.4885766651
12 edu_cat.maleSome Primary School (Grades 1-7) -0.4447284855
13 edu_cat.maleCompleted Primary School (Grade 7) -0.1868942957
14 edu_cat.maleSome Secondary School (Grades 8-10) -0.3185205153
15 edu_cat.maleCompleted Secondary School (Grade 10) -0.3999100980
16 edu_cat.maleCollege/Higher Education 0.1606673950
17 rcs(soc_sup_ps.female, 3)soc_sup_ps.female 0.0150021016
18 rcs(soc_sup_ps.female, 3)soc_sup_ps.female' -0.0115538098
19 rcs(soc_sup_ns.female, 3)soc_sup_ns.female -0.0370732576
20 rcs(soc_sup_ns.female, 3)soc_sup_ns.female' 0.0575217151
std.error statistic df p.value
1 1.3124172358 -0.30891420 595.2464 0.75749489
2 0.4006933936 -0.48508637 830.7336 0.62774297
3 0.0008159059 0.77376967 764.4603 0.43930623
4 0.0009997004 -0.46192410 665.5631 0.64428670
5 0.0412591653 0.44379338 624.4983 0.65734570
6 0.0602865086 0.01902215 537.8972 0.98483049
7 0.1587459455 1.35634812 146.8537 0.17707007
8 0.3051595733 2.32487375 363.4296 0.02062893
9 0.2772682392 2.18728433 532.9234 0.02915560
10 0.4914817981 1.59787495 715.4602 0.11051233
11 0.8568396612 0.57020781 480.6919 0.56880317
12 0.2214297901 -2.00844017 347.3996 0.04537008
13 0.2885136350 -0.64778323 521.6840 0.51741013
14 0.3175786969 -1.00296562 514.2414 0.31634915
15 0.3936528917 -1.01589524 371.1486 0.31034106
16 0.3765849591 0.42664315 581.6784 0.66979711
17 0.0373706931 0.40144028 366.4925 0.68832997
18 0.0584066884 -0.19781655 450.4450 0.84327783
19 0.0305405587 -1.21390240 344.0862 0.22561753
20 0.0437286979 1.31542255 426.4860 0.18907466# GEE with MD correction with geesmv package
prim_alt_aa_gee_check <- complete(imp_prim_aim2, "all") %>%
map(~ GEE.var.md(formula = ua_form,
family = binomial,
id = "clinic.female",
data = .x,
corstr = "exchangeable")) (Intercept) group.femaleIntervention
-0.5042467 0.1546620 (Intercept) group.femaleIntervention
-0.4693350 0.1567008 (Intercept) group.femaleIntervention
-0.50424665 0.07061067 (Intercept) group.femaleIntervention
-0.46065152 0.07387855 (Intercept) group.femaleIntervention
-0.4954913 0.1551655 (Intercept) group.femaleIntervention
-0.469334957 0.007357638 (Intercept) group.femaleIntervention
-0.4780358 0.1005837 (Intercept) group.femaleIntervention
-0.4260844 0.2594075 (Intercept) group.femaleIntervention
-0.44333538 -0.01864194 (Intercept) group.femaleIntervention
-0.4519851 0.1393510 (Intercept) group.femaleIntervention
-0.42608440 0.03931142 (Intercept) group.femaleIntervention
-0.40889564 -0.04361483 (Intercept) group.femaleIntervention
-0.4433354 0.1765965 (Intercept) group.femaleIntervention
-0.46065152 0.08319945 (Intercept) group.femaleIntervention
-0.4693350 0.1290092 (Intercept) group.femaleIntervention
-0.4780358 0.1746054 (Intercept) group.femaleIntervention
-0.4347019 0.2135806 (Intercept) group.femaleIntervention
-0.46933496 0.03569897 (Intercept) group.femaleIntervention
-0.48675444 0.02477712 (Intercept) group.femaleIntervention
-0.5218147 0.1443627 (Intercept) group.femaleIntervention
-0.4780358 0.1837963 (Intercept) group.femaleIntervention
-0.40889564 0.05003729 (Intercept) group.femaleIntervention
-0.5042467 0.1360993 (Intercept) group.femaleIntervention
-0.4088956 0.1421567 (Intercept) group.femaleIntervention
-0.4260844 0.1684896 (Intercept) group.femaleIntervention
-0.4347019 0.1588078 (Intercept) group.femaleIntervention
-0.45198512 0.08383778 (Intercept) group.femaleIntervention
-0.4867544 0.1833240 (Intercept) group.femaleIntervention
-0.4347019 0.2135806 (Intercept) group.femaleIntervention
-0.4433354 0.2131104 (Intercept) group.femaleIntervention
-0.47803580 0.03497228 (Intercept) group.femaleIntervention
-0.47803580 0.05380842 (Intercept) group.femaleIntervention
-0.4088956 0.1330015 (Intercept) group.femaleIntervention
-0.4693350 0.1290092 (Intercept) group.femaleIntervention
-0.4519851 0.1577457 (Intercept) group.femaleIntervention
-0.45198512 0.07453305 (Intercept) group.femaleIntervention
-0.460651524 0.008141052 (Intercept) group.femaleIntervention
-0.469335 0.302658 (Intercept) group.femaleIntervention
-0.4433354 0.1030096 (Intercept) group.femaleIntervention
-0.4780358 0.1561848 (Intercept) group.femaleIntervention
-0.4954913 0.1273439 (Intercept) group.femaleIntervention
-0.4174824 0.1232429 (Intercept) group.femaleIntervention
-0.5130210 0.1541627 (Intercept) group.femaleIntervention
-0.4519851 0.1760910 (Intercept) group.femaleIntervention
-0.4780358 0.1191775 (Intercept) group.femaleIntervention
-0.4693350 0.2845412 (Intercept) group.femaleIntervention
-0.4954913 0.3106975 (Intercept) group.femaleIntervention
-0.5130210 0.3009943 (Intercept) group.femaleIntervention
-0.4606515 0.1203257 (Intercept) group.femaleIntervention
-0.4433354 0.1490959 # GEE with my function, can test they are the same with the unadjusted formula
prim_alt_aa_geeglm_check <- complete(imp_prim_aim2, "all") %>%
map(~ geeglm.var.md(formula = ua_form,
family = binomial,
data = .x,
corstr = "exchangeable"))
# check that true to three decimal places
map2_df(prim_alt_aa_gee_check,
prim_alt_aa_geeglm_check,
~ round(.x$cov.beta, 3) == round(.y$cov.beta, 3))# A tibble: 2 × 50
`1` `2` `3` `4` `5` `6` `7` `8` `9` `10` `11` `12` `13`
<lgl> <lgl> <lgl> <lgl> <lgl> <lgl> <lgl> <lgl> <lgl> <lgl> <lgl> <lgl> <lgl>
1 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
2 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
# ℹ 37 more variables: `14` <lgl>, `15` <lgl>, `16` <lgl>, `17` <lgl>,
# `18` <lgl>, `19` <lgl>, `20` <lgl>, `21` <lgl>, `22` <lgl>, `23` <lgl>,
# `24` <lgl>, `25` <lgl>, `26` <lgl>, `27` <lgl>, `28` <lgl>, `29` <lgl>,
# `30` <lgl>, `31` <lgl>, `32` <lgl>, `33` <lgl>, `34` <lgl>, `35` <lgl>,
# `36` <lgl>, `37` <lgl>, `38` <lgl>, `39` <lgl>, `40` <lgl>, `41` <lgl>,
# `42` <lgl>, `43` <lgl>, `44` <lgl>, `45` <lgl>, `46` <lgl>, `47` <lgl>,
# `48` <lgl>, `49` <lgl>, `50` <lgl># now final version ## HIGHLY DEPENDED ON INCLUDED COVARIATES AND OUTCOME FOR SINGULARITY ##
## definitely have to remove who_stage, may have to remove district as well, but outcome dependent
## COULD NOT PUT DISTRICT BACK IN ##
# of note, ID is inherent to function (see above)
prim_alt_aa_md <- complete(imp_prim_aim2, "all") %>%
map(~ geeglm.var.md(formula = Y_imp ~ group.female + rcs(day.female, 3) + rcs(age.female, 3) +
+ edu_cat.female + edu_cat.male + rcs(soc_sup_ps.female, 3) + rcs(soc_sup_ns.female, 3),
family = binomial,
data = .x,
corstr = "exchangeable"))
# need to write function for GEEmd list $cov.beta[2] is what I'm interested in, since I only care about the intervention estimate :)
# using Rubin's rules from https://bookdown.org/mwheymans/bookmi/rubins-rules.html
# if m is the number of imputations
# vw = 1/m * sum(se^2)
# vb = sum(pooled beta - beta in each imputed dataset)^2 / (m - 1)
# vtotal = vw + vb+ vb/m
# se = sqrt(vtotal)
# m
m <- length(prim_alt_aa_md)
# first calculate vw
vw <- (1 / m) * sum(map_df(prim_alt_aa_md, ~ .x$cov.beta[2]))
# then calculate vb
num_vb <- map_df(prim_alt_aa, ~ .x$coefficients[2]) %>% # pull out intervention coefficients from imputed datasets
mutate(overall.coef = prim3$estimate[2],
diff = overall.coef - group.femaleIntervention) %>% # add difference between overal coefficient and coefficient from each dataset
pull(diff) %>% sum() # pull out that column and sum it
vb <- num_vb / (m - 1)
# now calculate vtotal
vtotal <- vw + vb + (vb / m)
se_prim3 <- sqrt(vtotal)Complete Case Adjusted Analysis
# unclear why this one was not working...removed WHO stage...then worked later...
prim_cc <- complete(imp_prim_aim2, "all") %>%
map(~ mm(mean.formula = cc_form,
lv.formula = ~ 1,
id = clinic.female,
q = 120,
data = subset(.x, !is.na(Y_raw)))) # now use my own pooling and summary functions to create output value
prim4 <- prim_cc %>%
pool.mm(., sum(!is.na(complete(imp_prim_aim2, 1)$Y_raw))) %>%
summary
prim4 term estimate
1 (Intercept) -1.4570347172
2 group.femaleIntervention 0.1507012647
3 rcs(day.female, 3)day.female 0.0005546062
4 rcs(day.female, 3)day.female' -0.0002639904
5 rcs(age.female, 3)age.female 0.0323793394
6 rcs(age.female, 3)age.female' -0.0234806337
7 district.femaleInhassunge 1.4173458800
8 district.femaleNamacurra 0.7004564181
9 district.femaleMocubela 1.7008674378
10 district.femaleMaganja da Costa 1.2274838244
11 district.femaleGilé 1.9408847705
12 district.femaleQuelimane -1.4838970700
13 who_stageII 0.1674695789
14 who_stageIII 0.1287877320
15 who_stageIV -4.8628311013
16 edu_cat.femaleSome Primary School (Grades 1-7) 0.4227559501
17 edu_cat.femaleCompleted Primary School (Grade 7) 0.8459876725
18 edu_cat.femaleSome Secondary School (Grades 8-10) 0.6846997608
19 edu_cat.femaleCompleted Secondary School (Grade 10) 0.6659271010
20 edu_cat.femaleCollege/Higher Education 0.2777566144
21 edu_cat.maleSome Primary School (Grades 1-7) -0.3796165408
22 edu_cat.maleCompleted Primary School (Grade 7) -0.3104841148
23 edu_cat.maleSome Secondary School (Grades 8-10) -0.3703053908
24 edu_cat.maleCompleted Secondary School (Grade 10) -0.5738272811
25 edu_cat.maleCollege/Higher Education 0.1139236450
26 rcs(soc_sup_ps.female, 3)soc_sup_ps.female 0.0249141115
27 rcs(soc_sup_ps.female, 3)soc_sup_ps.female' -0.0208620692
28 rcs(soc_sup_ns.female, 3)soc_sup_ns.female -0.0663701467
29 rcs(soc_sup_ns.female, 3)soc_sup_ns.female' 0.1018826951
std.error statistic df p.value
1 1.533788e+00 -0.949958252 740.3676 0.34244345
2 4.411362e-01 0.341620735 772.3038 0.73272926
3 6.921064e-04 0.801330844 754.2847 0.42319249
4 9.277556e-04 -0.284547398 758.1060 0.77606862
5 4.611193e-02 0.702190035 721.7399 0.48278705
6 6.346779e-02 -0.369961402 692.9871 0.71152429
7 6.286780e-01 2.254486171 772.5158 0.02444494
8 5.835483e-01 1.200340116 773.0151 0.23037499
9 6.724086e-01 2.529514663 773.6060 0.01161955
10 8.231657e-01 1.491174583 773.0828 0.13632357
11 7.545256e-01 2.572324649 773.8087 0.01028722
12 1.525528e+00 -0.972710388 774.6370 0.33100090
13 2.910009e-01 0.575495037 510.4183 0.56521016
14 4.991014e-01 0.258039219 285.3230 0.79656266
15 7.158214e+02 -0.006793358 774.9301 0.99458147
16 2.371266e-01 1.782827757 674.3124 0.07506412
17 3.702987e-01 2.284608964 771.8912 0.02260633
18 3.454609e-01 1.981989062 770.5920 0.04783534
19 5.065607e-01 1.314604832 768.5908 0.18903488
20 6.825285e-01 0.406952414 771.8981 0.68415576
21 3.056707e-01 -1.241913398 727.6638 0.21466864
22 3.594880e-01 -0.863684223 767.4916 0.38803121
23 3.483998e-01 -1.062874903 768.3787 0.28817261
24 4.258815e-01 -1.347387097 768.2262 0.17825285
25 4.539635e-01 0.250953322 766.4765 0.80191746
26 4.940939e-02 0.504238415 562.4570 0.61429131
27 7.207821e-02 -0.289436548 460.6584 0.77237755
28 4.433363e-02 -1.497061059 409.2449 0.13514820
29 5.582195e-02 1.825136957 378.8863 0.06876754Likely Case Missing Outcomes “No”
# unclear why this one was not working...removed District
prim_lc <- complete(imp_prim_aim2, "all") %>%
map(~ mm(mean.formula = lc_form,
lv.formula = ~ 1,
id = clinic.female,
q = 120,
data = .x)) # now use my own pooling and summary functions to create output value
prim5 <- prim_lc %>%
pool.mm(., nrow(complete(imp_prim_aim2, 1))) %>%
summary
prim5 term estimate
1 (Intercept) -0.8445215883
2 group.femaleIntervention -0.2459495318
3 rcs(day.female, 3)day.female 0.0007939286
4 rcs(day.female, 3)day.female' -0.0011392050
5 rcs(age.female, 3)age.female 0.0235606768
6 rcs(age.female, 3)age.female' -0.0133606863
7 who_stageII 0.1451291011
8 who_stageIII 0.0688484810
9 who_stageIV -4.6395030160
10 edu_cat.femaleSome Primary School (Grades 1-7) 0.3683299288
11 edu_cat.femaleCompleted Primary School (Grade 7) 0.7834160233
12 edu_cat.femaleSome Secondary School (Grades 8-10) 0.7489570836
13 edu_cat.femaleCompleted Secondary School (Grade 10) 0.6802095187
14 edu_cat.femaleCollege/Higher Education 0.1972430601
15 edu_cat.maleSome Primary School (Grades 1-7) -0.2889235420
16 edu_cat.maleCompleted Primary School (Grade 7) -0.1354598703
17 edu_cat.maleSome Secondary School (Grades 8-10) -0.3122756309
18 edu_cat.maleCompleted Secondary School (Grade 10) -0.3659986054
19 edu_cat.maleCollege/Higher Education 0.1527163922
20 rcs(soc_sup_ps.female, 3)soc_sup_ps.female 0.0261346981
21 rcs(soc_sup_ps.female, 3)soc_sup_ps.female' -0.0367424381
22 rcs(soc_sup_ns.female, 3)soc_sup_ns.female -0.0508857589
23 rcs(soc_sup_ns.female, 3)soc_sup_ns.female' 0.0753121948
std.error statistic df p.value
1 1.343864e+00 -0.628427799 895.7609 0.52988393
2 4.205109e-01 -0.584882606 927.5040 0.55876890
3 6.095411e-04 1.302502205 919.6686 0.19307084
4 8.404228e-04 -1.355514150 921.7894 0.17558583
5 4.356506e-02 0.540815944 406.3103 0.58893081
6 6.764819e-02 -0.197502493 373.8765 0.84354166
7 2.706694e-01 0.536185865 540.5568 0.59205067
8 4.673893e-01 0.147304363 308.4207 0.88298804
9 7.437093e+02 -0.006238329 928.9133 0.99502391
10 2.168751e-01 1.698350325 846.1963 0.08980928
11 3.335068e-01 2.349025390 923.4801 0.01903253
12 3.191643e-01 2.346619424 919.8250 0.01915594
13 4.565375e-01 1.489931216 924.5109 0.13658339
14 5.762939e-01 0.342261243 927.5233 0.73223184
15 2.770999e-01 -1.042669133 880.2414 0.29738783
16 3.249328e-01 -0.416885810 926.5094 0.67685850
17 3.133119e-01 -0.996692641 923.1754 0.31917489
18 3.751799e-01 -0.975528224 923.8117 0.32955368
19 4.016749e-01 0.380198952 923.5275 0.70388516
20 5.215794e-02 0.501068479 234.6279 0.61679269
21 6.926232e-02 -0.530482368 206.3013 0.59634808
22 4.346095e-02 -1.170838501 296.4773 0.24260390
23 5.346939e-02 1.408510389 265.6685 0.16014921Worst-case scenario
# took out district
prim_wc <- complete(imp_prim_aim2, "all") %>%
map(~ mm(mean.formula = wc_form,
lv.formula = ~ 1,
id = clinic.female,
q = 120,
data = .x)) # now use my own pooling and summary functions to create output value
prim6 <- prim_wc %>%
pool.mm(., nrow(complete(imp_prim_aim2, 1))) %>%
summary
prim6 term estimate
1 (Intercept) -0.6828212762
2 group.femaleIntervention -0.6677408608
3 rcs(day.female, 3)day.female 0.0006168625
4 rcs(day.female, 3)day.female' -0.0002752680
5 rcs(age.female, 3)age.female 0.0157004069
6 rcs(age.female, 3)age.female' 0.0072657606
7 who_stageII 0.0554077041
8 who_stageIII -0.1157159897
9 who_stageIV -4.6725738537
10 edu_cat.femaleSome Primary School (Grades 1-7) 0.2655548859
11 edu_cat.femaleCompleted Primary School (Grade 7) 0.7742952169
12 edu_cat.femaleSome Secondary School (Grades 8-10) 0.5902249688
13 edu_cat.femaleCompleted Secondary School (Grade 10) 0.3071245658
14 edu_cat.femaleCollege/Higher Education 0.7198684126
15 edu_cat.maleSome Primary School (Grades 1-7) -0.0730508025
16 edu_cat.maleCompleted Primary School (Grade 7) 0.0824390290
17 edu_cat.maleSome Secondary School (Grades 8-10) -0.0400800684
18 edu_cat.maleCompleted Secondary School (Grade 10) 0.0366083106
19 edu_cat.maleCollege/Higher Education 0.2864409796
20 rcs(soc_sup_ps.female, 3)soc_sup_ps.female 0.0156598602
21 rcs(soc_sup_ps.female, 3)soc_sup_ps.female' -0.0036320302
22 rcs(soc_sup_ns.female, 3)soc_sup_ns.female -0.0361347281
23 rcs(soc_sup_ns.female, 3)soc_sup_ns.female' 0.0467407826
std.error statistic df p.value
1 1.299324e+00 -0.525520540 878.8377 0.59935394
2 4.254296e-01 -1.569568334 926.2362 0.11685710
3 6.041015e-04 1.021123932 814.5721 0.30749889
4 8.180280e-04 -0.336501941 831.2909 0.73657731
5 3.998337e-02 0.392673420 549.3111 0.69471303
6 6.058980e-02 0.119917221 593.7301 0.90458929
7 2.583423e-01 0.214473969 631.5303 0.83024670
8 4.570063e-01 -0.253204369 322.8386 0.80027145
9 7.431200e+02 -0.006287779 928.9133 0.99498446
10 2.112708e-01 1.256940680 700.6636 0.20919396
11 3.258375e-01 2.376323236 921.9984 0.01768952
12 3.084880e-01 1.913283096 917.7892 0.05602270
13 4.598560e-01 0.667871160 922.9217 0.50438281
14 5.547118e-01 1.297734141 927.4585 0.19470137
15 2.781582e-01 -0.262623206 672.2671 0.79292147
16 3.122689e-01 0.264000086 924.9538 0.79183866
17 3.004675e-01 -0.133392347 920.5315 0.89391225
18 3.600244e-01 0.101682848 922.4811 0.91903051
19 3.949172e-01 0.725318995 923.8326 0.46844004
20 4.346453e-02 0.360290545 525.0941 0.71877475
21 5.737313e-02 -0.063305426 491.9568 0.94954901
22 3.914165e-02 -0.923178521 435.1839 0.35642574
23 4.537044e-02 1.030203432 596.5561 0.30333201COVID-19 Sensitivity Analyses
Load(imp_prim_aim2_cov)
# first need to make clinic.female a factor in all iterations
long <- complete(imp_prim_aim2_cov, action = "long", include = TRUE)
long$clinic.female <- factor(long$clinic.female)
# now make it back into a mids object
imp_prim_aim2_cov <- as.mids(long)Unadjusted analysis
cov_ua <- complete(imp_prim_aim2_cov, "all") %>%
map(~ mm(mean.formula = ua_form,
lv.formula = ~ 1,
id = clinic.female,
q = 120,
data = .x))
# nw use my own pooling and summary functions to create output value
cov1 <- cov_ua %>%
pool.mm(., nrow(complete(imp_prim_aim2_cov, 1))) %>%
summary
cov1 term estimate std.error statistic df p.value
1 (Intercept) -0.3469851 0.2553812 -1.358695 528.5033 0.1748230
2 group.femaleIntervention -0.4217567 0.3805475 -1.108289 500.1909 0.2682696Adjusted analysis
cov_aa <- complete(imp_prim_aim2_cov, "all") %>%
map(~ mm(mean.formula = a_form,
lv.formula = ~ 1,
id = clinic.female,
q = 120,
data = .x)) # nw use my own pooling and summary functions to create output value
cov2 <- cov_aa %>%
pool.mm(., nrow(complete(imp_prim_aim2_cov, 1))) %>%
summary
cov2 term estimate
1 (Intercept) -4.041362447
2 group.femaleIntervention -0.132485930
3 rcs(day.female, 3)day.female 0.001706540
4 rcs(day.female, 3)day.female' -0.001841392
5 rcs(age.female, 3)age.female 0.090715943
6 rcs(age.female, 3)age.female' -0.080303587
7 district.femaleInhassunge 0.966838698
8 district.femaleNamacurra 0.561218881
9 district.femaleMocubela 1.501480215
10 district.femaleMaganja da Costa 0.028024492
11 district.femaleGilé 1.826525584
12 district.femaleQuelimane -1.218430829
13 who_stageII 0.153380443
14 who_stageIII 0.123603422
15 who_stageIV -6.764927305
16 edu_cat.femaleSome Primary School (Grades 1-7) 0.127091808
17 edu_cat.femaleCompleted Primary School (Grade 7) 1.027787228
18 edu_cat.femaleSome Secondary School (Grades 8-10) 0.611753680
19 edu_cat.femaleCompleted Secondary School (Grade 10) 1.468748423
20 edu_cat.femaleCollege/Higher Education -0.936038150
21 edu_cat.maleSome Primary School (Grades 1-7) -0.473423618
22 edu_cat.maleCompleted Primary School (Grade 7) 0.094254171
23 edu_cat.maleSome Secondary School (Grades 8-10) -0.152571110
24 edu_cat.maleCompleted Secondary School (Grade 10) -0.237212311
25 edu_cat.maleCollege/Higher Education 0.728794744
26 rcs(soc_sup_ps.female, 3)soc_sup_ps.female 0.059818832
27 rcs(soc_sup_ps.female, 3)soc_sup_ps.female' 0.003769539
28 rcs(soc_sup_ns.female, 3)soc_sup_ns.female -0.036803354
29 rcs(soc_sup_ns.female, 3)soc_sup_ns.female' 0.063377778
std.error statistic df p.value
1 2.259612e+00 -1.78852036 362.7967 0.074526374
2 3.687148e-01 -0.35931822 466.5858 0.719519604
3 2.032912e-03 0.83945592 353.2604 0.401781238
4 2.733595e-03 -0.67361544 324.1859 0.501035798
5 7.174791e-02 1.26437053 302.3053 0.207071004
6 8.294800e-02 -0.96811960 299.6950 0.333765260
7 5.004743e-01 1.93184471 515.9041 0.053925947
8 4.762133e-01 1.17850309 522.1672 0.239132861
9 5.800097e-01 2.58871575 340.8404 0.010046080
10 7.257608e-01 0.03861395 301.8155 0.969223708
11 5.877786e-01 3.10750610 530.9507 0.001987725
12 1.286651e+00 -0.94697853 534.0832 0.344077973
13 3.984800e-01 0.38491382 287.8261 0.700585535
14 6.863191e-01 0.18009615 147.9590 0.857323513
15 8.547888e+02 -0.00791415 539.9569 0.993688411
16 3.271503e-01 0.38848142 362.6301 0.697887777
17 4.873251e-01 2.10903813 412.2815 0.035543756
18 4.655255e-01 1.31411417 454.8516 0.189470114
19 7.383326e-01 1.98927742 370.5628 0.047405398
20 9.909332e-01 -0.94460270 338.0331 0.345536416
21 4.500032e-01 -1.05204510 318.5097 0.293576346
22 5.039675e-01 0.18702430 363.6516 0.851745877
23 5.016844e-01 -0.30411772 357.2282 0.761215274
24 6.211124e-01 -0.38191527 304.5230 0.702790584
25 6.230094e-01 1.16979739 373.4668 0.242828512
26 8.493405e-02 0.70429745 154.9450 0.482304739
27 1.030744e-01 0.03657107 152.3509 0.970874889
28 6.944555e-02 -0.52995987 199.5187 0.596729011
29 8.347705e-02 0.75922394 176.0923 0.448733722